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ABSTRACT 

Astrophysical discs are warped whenever a misalignment is present in the system, or 
when a flat disc is made unstable by external forces. The evolution of the shape and 
mass distribution of a warped disc is driven not only by external influences but also by 
an internal torque, which transports angular momentum through the disc. This torque 
depends on internal flows driven by the oscillating pressure gradient associated with 
the warp, and on physical processes operating on smaller scales, which may include in- 
stability and turbulence. We introduce a local model for the detailed study of warped 
discs. Starting from the shearing sheet of Goldreich & Lynden-Bell, we impose the os- 
cillating geometry of the orbital plane by means of a coordinate transformation. This 
warped shearing sheet (or box) is suitable for analytical and computational treatments 
of fluid dynamics, magnetohydrodynamics, etc., and it can be used to compute the 
internal torque that drives the large-scale evolution of the disc. The simplest hydrody- 
namic states in the local model are horizontally uniform laminar flows that oscillate 
at the orbital frequency. These correspond to the nonlinear solutions for warped discs 
found in previous work by Ogilvie, and we present an alternative derivation and gen- 
eralization of that theory. In a companion paper we show that these laminar flows are 
often linearly unstable, especially if the disc is nearly Keplerian and of low viscosity. 
The local model can be used in future work to determine the nonlinear outcome of 
the hydrodynamic instability of warped discs, and its interaction with others such as 
the magnetorotational instability. 
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1 INTRODUCTION 

1.1 Astrophysical motivation 

Warped discs, in which the orbital plane varies with ra- 
dius, have many applications in astrophysics. They occur 
whenever a misalignment is present in the system, as in 
the classic scenario of a black hole whose spin axis does 
not coincide with the or bital axis of gas that is su pplied 
through an accretion disc ([Bardeen fc Pettersonlll975l ) . Vari- 
ants of this problem occur when the central object is a 
magnetized star or a close binary, interacting with the disc 
through magnetic or gravitational torques. A disc may also 
be warped by a companion object on an inclined orbit 
l|Papaloizou fc Terauemlirggsl ) ; this situation is found in suf- 
ficiently wide young binary stars and can occur in proto- 
planetary systems if a mutual inclination of the planet and 
disc is excited by many-body interactions. Even in systems 
that are initially coplanar, warps may arise spontaneously 
through the growth of ins tabilities , such as those involvin g 
tidal forces l|Lubo w' 1992), winds (S chandl fc M eyer '1994). 
radiation forces (fPri ngle 1996) or magnetic fields (Lai 1999i). 
Early studies of warped discs were motivated not only 



by theoretical proble ms such as misaligned accre tion on to 
a spinning black hole l|Bardeen fc Petterson|[l975l) . but also 
by observational discoveries such as the low-mass X-ray bi- 
nary Her X-1 (HZ Her), where the existence of a precess- 
ing disc tilte d out of th e binary plane was deduced from 
light curves (|Kat j 0*9731 ) . Much later, observations of wa- 
ter masers in the galaxy NGC 4258 (M 1 06) revealed a 
warp ed disc around the central black hole (iMivoshi et al.l 
ll995l V There are by no w multiple examples of X-ray binaries 
(e.g. Kotze fc Charlej|2012l ) and active galactic nuclei (e.g. 
lGreenhil]||2005l ) that may have similar properties to these 
systems. Recent interest in warped discs has focused mainly 
on applications to accreting b l ack h o les and to protoplan- 
etary syste ms. iNixon fc Kind (|2012D . H ixon. King fc Pricg 
(j2012h and lNixon et al.l l|2012h have argued that, in signifi- 
cantly misaligned accretion on to a spinning black hole, the 
disc breaks into rings that can precess independently, and 
the accretion rate is greatly enhanced, jjoucart fc Lai (201lj) 
have calculated the warping of a protoplanetary disc that is 
tilted with respect to the spin axis of a magnetized central 
star, and have investigated the consequences of this dynam- 
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ics for tfie spin-orbit misalignment of extrasolar planetary 
systems. 

1.2 Theoretical and computational background 

Fundamental theoretical studies of warped discs have mostly 
aimed to derive equations that govern the evolution of the 
shape of the disc and, in some cases, its mass distribution. 
There is also an extensive literature that applies the theory 
of warped discs, often in a simplified form, to astrophysical 
systems. 

Early versi ons of evolutionary equations for the shape of 
a warped disc |Bardeen fc Petters"oiil 19751: |Pettersonlll977l . 
ll978l : lHatchett Begehnan, fc Sarazinl ligSl*) differed slightly 
from each other but all suggested that the warp would diffuse 
on a viscous timescale. They all turned out to be incorrect, 
mainly because the internal flows driven by the oscillating 
press ure gradient in a warped d isc had not been consid- 
ered. IPagaloizou^&^^ringla (1983) provided the first consis- 
te nt linear theory for visc ous Keplerian discs (summarized 
bv lKumar fc Pringldl 19851 ). and found that the warp diffuses 
on a timescale that is shorter than the viscous timescale by a 
factor of order a^, where a <C 1 is the Shakura-Sunyaev vis- 
cosity parameter. (In this context, viscosity is usually taken 
to represent the effects of unresolved phys ical processes such 
as sm all-scale turbulence.) Subsequently, IPapaloizou fc LinI 
l|l995 h showed that a transition from diffusive to wavelike 
propagation occurs in a Keplerian disc when a is less than 
the angular s emith ickness H/r of the disc. 

lOgilvid (|l999l ) derived a fully nonlinear theory for 
the diffusive regime in Keplerian discs and also for 
non-Keplerian discs. While the resulting evolutionary 
equations are simila r in form to tho s e su ggested by 
IPaoaloizou fc Pringld (|l983l ) and IPringld (|l992l ). the anal- 
ysis also provides a means to calculate the torque coeffi- 
cients in those equations as functions of the amplitude of the 
warp and other relevant parameters. The thermal physics of 
w arped discs in t he nonlinear reg ime was taken into account 
bv lOgilvid (|2000D . More recent Iv. lOgiM^ (|2006l ) showed how 
the wavelike regime for Keplerian discs is modified by weak 
nonlinearity and dispersion: solitary bending waves would 
be possible if the adiabatic exponent 7 were to exceed 3, 
but otherwise the dominant weakly nonlinear effect is to 
enhance the dispersive spreading of a bending wave. 

Global numerical simulations of warped discs are very 
demanding because of the ranges of length-scales and time- 
scales that are involved in a thin disc. While there have 
been some grid-based simulations of warped discs, the ma- 
jority of studies have used smoothed particle hydrody- 
namics (SPH), which is well suited to the complicated 
geometry, although less good for resolving small scales. 
The evolution of simple warps in SPH simulations was 
measured and compared w it h theoretical expe c tation s by 
Nelson fc Papaloizoul (1 19991 ) . iLodato fc Pringld (I2007D and 
Lodato fc Pried (|2010^ . The last paper, in p articul a r, con - 
firms some aspects of the nonlinea r t heory ofiOgilvid (1 19991). 
Previo usly, iLarwood et al.l l| 19961 ) , iLarwood fc Papaloizoul 
l|l997h and iLarwoodT i 19971 ) had studied the interaction of 
discs with bi nary companions on inclined orbit s; the re- 
cent work of IXiang-Gruess fc Papal oizou' (2013') involves 
planets on inclined orbits. iNelson fc P apaloizou (2000) in- 
cluded a post-Newtonian force within SPH to simulate 



tilted discs around spinning black holes. The tilting and 
warping of discs in binary stars by magnetic an d radi - 



atio n forces has been simulate d bv [ Murray et al.l (|2002l ) 
and iFoulkes. HasweU fc Murravl l|2006l . I2OI0I ). Simulations 



using grid-based methods have been applied to tilted 
discs around spinn ing black holes (iFragile fc Ann inos"2005l: 
[Fragile et al. I l2007l . [2009. ). where they reveal a complicated 
beh aviour of the accret i on str eams close to the event hori- 
zon. iFragner fc NelsonI (|2010l ) used rotating grids to com- 
pute precessing circumstellar discs in binary stars, obtaining 
general agreement with theoretical expectations, but finding 
that in extreme cases the disc tends to break into indepen- 
dently precessing rings. This type of behaviour has also been 
emphasized in the previously mentioned work by Nixon (e.g. 
iNixon et al.„20ia ). which uses SPH. 



1.3 Plan of this paper 

The main purpose of this paper is to introduce a local 
model for the detailed study of warped discs, including in- 
stability and turbulence. We first discuss the large-scale ge- 
ometry of a warped disc and show how the evolution of 
its shape and mass distribution is driven by an internal 
torque. We then use a circular reference orbit to construct 
a standard local model, equiv alent to the shearing sheet 
l|Goldreich fc Lvnden-Belilll965l ) and including vertical grav- 
ity. To take into account the oscillating local geometry of the 
orbital plane we then introduce a transformation to warped 
shearing coordinates and formulate the hydrodynamic equa- 
tions in this new system. A single dimensionless parameter 
defines the local properties of the warp in this model, and 
we show how to compute the internal torque that governs 
the large-scale evolution of the disc. The simplest hydro- 
dynamic states in the warped shearing box are horizontally 
uniform laminar flows that oscillate at the orbital frequency. 
We explore the properties of these laminar flows and the re- 
lated torques, which corre spond to the n onlinear solutions 
for wa rped discs found bv [ Ogilvid (Il999l ). In a companion 
paper (|Ogilvie fc LatteilbOia hereafter Paper II) we use the 
local model to analyse the linear hydrodynamic stability of 
the laminar flows and find widespread instability, which re- 
quires further investigation in future work. 



2 LARGE-SCALE GEOMETRY AND 
DYNAMICS OF A WARPED DISC 

In a thin astrophysical disc, the orbital motion is hyper- 
sonic and fluid elements follow ballistic trajectories to a first 
approximation. Around a spherical central mass, these tra- 
jectories are Keplerian orbits, which can have eccentricity 
and inclination. A general Keplerian disc involves smoothly 
nested orbits of variable eccentricity and inclination: it is 
both elliptical and warped. 

We consider a spherically symmetric gravitational po- 
tential in which circular orbital motion is possible in any 
plane containing the centre of the potential^ The dominant 
motion in a warped disc is orbital motion in a plane that 



^ Any small non-spherical component of the potential can be con- 
sidered to provide an external torque on the disc. 
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Figure 2. Edge-on view of discs with untwisted warps of constant 
amplitudes \ip\ = 0.01 (top), 0.1, 0.2, 0.5 and 1 (bottom). Each 
curve is a logarithmic spiral. 



Figure 1. Top: Example of an untwisted warped disc, viewed 
as a collection of tilted rings. As discussed in Section |3] a refer- 
ence orbit (red circle) is selected and used to construct a local 
Cartesian model of the disc. The local frame (blue box) is cen- 
tred on a point that follows the reference orbit, and therefore 
experiences a geometry that oscillates at the orbital frequency as 
the local orientation of the midplane of the disc tilts back and 
forth. The illustrated box corresponds to the local frame at or- 
bital phase according to our definitions, and its axes are aligned 
with the radius-dependent basis vectors {I, m, n} at this phase 
only. In this example the warp amplitude at the reference radius 
is \ip\ = 0.5. Bottom: Edge-on view of the warped disc, showing 
the local basis vectors. 



varies continuously with radius r and possibly with time t. 
In fact the motion need not be circular, but we will not con- 
sider eccentric warped discs in this paper. We assume that 
the time-dependence is slow so that the shape of the warp 
can be regarded as fixed on the orbital time-scale. There- 
fore the warped disc can be considered as a continuum of 
tilted rings (Fig. [T] top). Since astrophysical discs are of 
non-zero thickness, this structure can be considered to de- 
fine the warped midplane of the disc. 

The orbital angular velocity of the warped disc is 
f2(r, t) = fl{r) l{r, t), where the slowly evolving unit tilt vec- 
tor l{r,t) is everywhere normal to the local orbital plane. 
The rate of orbital shear is 



dft dfi , ^ dl 

: r— = r— I + ril — 
or dr or 



(1) 



where q — — dlnS7/dlnr is the usual dimensionless rate 
of orbital shear, equal to | for Keplerian orbits, \ip\ — 



\dl/dlnr\ is t he dimensionless warp amplitude defined by 
lOgilvidll 19991) . and m is a unit vector parallel to dl / dr and 
therefore orthogonal to L A right-handed orthonormal triad 
{I, m, n}, dependent on r and t but not on the orbital phase, 
is completed hy n = I x m (Fig. [TJ bottom) . Note that none 
of these vectors is normal to the disc, in general. 

The simplest type of warp is an 'untwisted' warp, by 
which we mean that, as in Fig. [T] the variation of the vec- 
tors I and m is confined to a plane, while the vector n is 
perpendicular to that plane and independent of r. In the 
language of celestial mechanics, the longitude of the ascend- 
ing node is independent of the semimajor axis. Although for 
clarity we choose untwisted warps for the purposes of illus- 
tration, our analysis is valid for general, twisted warps. In- 
deed, the parameter \^\ measures the local amplitude of the 
warp whether or not it is twisted. A local measure of twist 
would involve the second radial derivative of I, for exam- 
ple the triple scalar product of Z, dl/dlnr and d^l/d{hir)^ . 
(Some authors, however, have used the term 'twisted disc' 
as synonymous with 'warped disc'.) 

The significance of the parameter is illustrated in 
Fig. [21 where we present an edge-on view of discs with un- 
twisted warps corresponding to difi'erent values of \tp\. We 
will see in this paper and its companion that a warp of am- 
plitude \tp\ = 0.01, which might not be directly observable 
even with high-resolution imaging, can nevertheless have im- 
portant dynamical consequences. 

The large-scale dynamics of a warped disc is gov- 
erned by the conservation o f mass and an gular momentum 
jPapaloizou fc Pringle|[l983l : |Pringle|[l993 ). The usual aim 
of a theory of warped discs is to obtain a system of par- 
tial differential equations in r and t that govern the shape 
and mass distribution of the disc. For a thin disc in which 
the angular momentum is dominated by orbital motion, the 
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conservation laws in the absence of external influences can 
be written in the one- dimensional form 

dt dr 



0, 



|(M/^) + |:(p/^ + g) = o, 



(2) 
(3) 



where A4 (r, t) is the one-dimensional mass density of the 
disc (related to the surface density E through A4 — 27rrE 
and deflned such that the mass contained between any two 
radii is given by an integral J Mdr between those radii) , 
J-{r, t) is the outward radial flux of mass and h{r, t) = 
r^f2(r, = h{r)l{r,t) is the specific angular momentum. 
The outward radial fiux of angular momentum consists of 
two parts: J-{r,t)h{r,t) due to advection and Q{r,t) being 
an internal torque associated with internal flows and stresses 
due to viscosity, magnetic fields, turbulence, self-gravitation, 
etc. Subtracting h times equation ((5)) from equation we 
obtain the equation governing the shape of the disc. 



Mh—+T— + ^ = 0. 

at or or 

A scalar product with the unit vector I gives 



dr dr 



= 0, 



(4) 



(5) 



which determines F in terms of Q. Therefore our task is 
to determine the internal torque Q. It is the transport of 
angular momentum that drives the evolution of both the 
shape of the disc and its mass distribution. The reason for 
this is of course that, within the family of circular orbits, 
the specific angular momentum vector determines both the 
orbital plane and the orbital radius. 

This analysis is easily extended to allow for external 
forces by including a source term T, the external torque per 
unit radius, on the right-hand side of equation 

Although the internal torque Q at any radius in a thin 
disc involves an integral with respect to the azimuthal angle 
and is in this sense a global or large-scale quantity, we will 
see below that this integral naturally emerges in the form 
of a time-average in a local model that follows the orbital 
motion. 



3 LOCAL MODEL OF A WARPED DISC 

3.1 Geometry and particle dynamics 

A local model can be constructed around any reference point 
that is situated on the warped midplane of the disc and 
moves in a circular orbit of radius ro, where the orbital an- 
gular velocity is f2o = f2('"o)- A Cartesian coordinate system 
is set up with its origin at the moving reference point and 
with axes that rotate with the orbital motion, so that the x, 
y and z directions are always radial, azimuthal and vertical. 
This is the standard construction of a local model, as used 
in the shearing sheet or shearing box. 

The effective potential in this frame is the sum of the 
gravitational potential and the centrifugal potential arising 
from the rotation of the frame. When the effective poten- 
tial is expanded in a Taylor series about the reference point 
and the unimportant constant term is neglected, the domi- 




Figure 3. Time-dependent relation between the rotating hori- 
zontal basis vectors ex and ey of the local model and the non- 
rotating basis vectors mp and ng associated with the geometry 
of the warp at the reference radius. The vertical basis vectors 
and lo both point out of the page. 



nant terms are —qQoX^ + IQ.qz'^. The equations of particle 



dynamics in this model are therefore 
y + 2nox = 0, 



(6) 
(7) 
(8) 



which reduce to Hill's equations (without a satellite) in the 
Keplerian case q — ^. 

As stated above, we assume that the geometry of the 
warped disc can be regarded as fixed on the orbital time- 
scale relevant to the local model, although we discuss this as- 
sumption further in Section [4.41 below. The basis vectors as- 
sociated with the warp at the reference radius, {Iq, mo, no}, 
are therefore non-rotating, while the basis {ex,ey,ez} ro- 
tates as the reference point traverses its orbit. Without loss 
of generality, we choose the orbital phase of the reference 
point, or the origin of time, such that mo is in the radial 
direction at t = (as in Fig. [ij . Then the basis vectors are 
related by 



lo 



mo + ino = (e^, + iey)e' 



(9) 



(Fig. |3]). While the rate of orbital shear is stationary in a 
non-rotating frame, in the local frame it rotates in a negative 
sense according to 



So^t) = —qilo ez + \rlj\Qo[ex cos(f2ot) - 



i(f^ot)], (10) 



which follows from equations ([T} and In Fig. [T] the box 
can be thought of as following the red reference orbit, expe- 
riencing a local geometry that oscillates at the orbital fre- 
quency as the local orientation of the midplane of the disc 
tilts back and forth. 

The local representation of circular orbital motion in 
the xy plane is the linear motion y — —qflox with x — 
constant. To this can be added a free vertical oscillation 
z — 'Re{Z e~^^°'^). The amplitude of the complex number 
Z is related to the (small) inclination of the orbit with re- 
spect to the xy plane, while the phase of Z is related to 
the longitude of the ascending node ; the c omplex tilt vari- 
able (used by iPapaloizou fc Pringle|[l98^ and many other 
authors), measured with respect to the xy plane, would be 
W = -Z/ro. 
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An orbital motion with a smooth warp can be repre- 
sented locally by allowing Z to vary linearly with x. Since 
the reference point is on the warped midplane of the disc, Z 
vanishes at x = 0. With the choice of orbital phase explained 
above, we have Z = — |-i/)|a;, and so z = —\tp\ coa{Qot)x. This 
motion corresponds to the velocity field 



X — u = —qQ,QX Sy + |i/;|r2o sin(r2oi)2; 



(11) 



Thus, while a particle at the origin remains there (corre- 
sponding to the reference orbit), a particle with a positive 
value of X lags behind and also oscillates up and down at the 
orbital frequency (corresponding to an orbit slightly larger 
than the reference orbit). 

Later we will need to calculate the internal torque in the 
local model. To prepare the way for this, we consider here 
the specific angular momentum, which in the local model 
can be regarded as 



(12) 



The first two vectors involve non-radial motions combined 
with the long radial lever arm, while the third vector in- 
volves the large azimuthal motion combined with a vertical 
lever arm. (The 2Qox term represents the contribution to 
the azimuthal motion from the rotation of the frame of ref- 
erence.) If we had also included the large azimuthal motion 
combined with the long radial lever arm, we would have 
obtained an additional term rgQa e.z ~ ho, which is con- 
stant and large compared to the terms considered here. For 
a general particle motion governed by equations ((6])-((8|, it 
can be verified the vector (|12p is constant in a non-rotating 
frame. Its components in the local frame, however, are not 
constant. 

For the orbital motion (jlip associated with a warped 
disc, the specific angular momentum evaluates to 



ro(So + 2no)x= ( ^ 



(13) 



and therefore agrees with a linear local approximation to the 
orbital angular momentum. Again, if we had also included 
the large azimuthal motion combined with the long radial 
lever arm, we would have obtained an additional term ho, 
which is constant and large compared to the terms consid- 
ered here. 



3.2 Fluid dynamics 

So far we have considered the motion of test particles. The 
local hydrodynamic solutions for a warped disc are more 
complicated than the particle motion (|lip because of pres- 
sure gradients. These solutions generally require numerical 
calculations and are most easily obtained by introducing a 
coordinate transformation that accounts for the warped or- 
bital motion. Before doing so we establish the equations to 
be solved. 

The basic equations for an ideal fluid in the local model 
(neglecting magnetic fields and self-gravity) are 



2 1 

T)Ux — 2^loUy = 2qfloX dxP, 



Duy + 2floUx = —-dyp, 



(14) 
(15) 



Duz — —QIz — -dzP, 



Dp = —p{dxUx -f dyUy -I- dzUz), 
where 

D = 9t -I- Uxdx + Uydy + Uzdz 



(16) 
(17) 

(18) 



is the Lagrangian derivative, u is the velocity, p is the den- 
sity and p is the pressure. These are the standard equa- 
tions for hydrodynamics in the s hearing sheet or box (e.g. 
lHawlev. Gammie fc Balbus|[l995h . For simplicity, we con- 
sider an isothermal gas for which 



(19) 



where Cs = constant is the isothermal sound speed. In terms 
of the pseudo-enthalpy 

h — c^lnp + constant, 



we then have 
Dux — 2floUy 
Duy + 2QoUx 



2qD,oX 



dxh, 



Duz = —D.oz — dzh, 

Dh = —c^idxUx + dyUy -f dzUz 



(20) 

(21) 
(22) 
(23) 
(24) 



The alternative equations for adiabatic flow are developed 
in Appendix 1X1 

Note that our local model correctly includes the vertical 
gravitational acceleration —Qq^ deriving from a spherically 
symmetric potential. Without this term the model would be 
incapable of representing a warped disc and the dynamics 
described in this paper would not occur. 

If the disc is unwarped, the orbital plane everywhere 
coincides with the xy plane and the local representation of 
circular orbital motion is u = —qfloxsy. Equations (I21|| - 
(|24|l are then satisfied when hydrostatic equilibrium, dzh = 



— Jlg^, is also imposed. 

If the disc is warped at the reference radius, then we 
might expect the velocity field (|lip corresponding to the 
warped orbital motion to satisfy the hydrodynamic equa- 
tions (|21|| - (|24|) . together with some version of hydrostatic 
equilibrium. However, this is not the case; the problem lies 
with equation (|23|l . If hydrostatic equilibrium, dzh = —QqZ, 
is imposed, then there is an unbalanced, a;-dependent accel- 
eration on the left-hand side of equation (|23p . If the pressure 
is allowed to vary with x to compensate for this term, then 
equation (|21|l is disrupted. The actual hydrodynamic solu- 
tions must be more complicated than equation (|ll|l in order 
to account for the interplay of these additional forces. The 
simplest (laminar) versions of these solutions are derived in 
Section [4] below. 



3.3 Warped shearing coordinates 

We now adapt the local model to incorporate the oscillating 
local geometry of the warp explicitly. This will allow us to 
find the simplest hydrodynamic states and to formulate a 
theoretical and computational model for the further study 
of warped discs. 

We introduce new coordinates 



t = t, 



(25) 
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IY'I cos t 



Figure 4. Illustration of warped shearing coordinates. The dotted blue grid represents the Cartesian coordinates of the standard local 
approximation or shearing box, while the solid black grid represents the warped shearing coordinates. The shear in the xy plane (left) 
increases linearly with time, although in a computational model the grid should be remapped periodically by resetting the origin of time. 
The shear in the xz plane (right) oscillates at the orbital frequency and is proportional to the amplitude of the warp. 



y' = y + qflotx, 



(26) 
(27) 
(28) 



The hydrodynamic equations are therefore transformed 



that foUow the warped orbital motion: a particle on an in- 
clined circular orbit would have constant x' , y' and z' . We 
define the orbital phase 

T = Hot' = not. (29) 
Partial derivatives transform according to 

dt — d't + qi^oxd'y — liplflosinrxd'^, (30) 

dx = di,. + qr d'y + l-tplcosT d'z, (31) 

dy = d'y, (32) 

d.^d'., (33) 
so that 

T> = d't +Vxd'^ + {vy +qTVx)d'y + (uz + |^/)| cosrwa;)^^, (34) 
where 

v^=u^, (35) 

Vy — Uy + q^lox, (36) 

Vz ~ Uz — Itpl^o siiiT X (37) 

are the relative velocity components. Thus, if v = (which 
is not a solution of the equations), the fluid follows the pre- 
scribed warped orbital motion and the velocity u, which cor- 
responds to the rate of change of the Cartesian coordinates 
{x, y, z), is non-zero. (We should not call u the absolute ve- 
locity, however, because it is measured in a rotating frame 
of reference.) We continue to refer vector components to the 
Cartesian basis {e^, By, Bz}. The grid of warped shearing co- 
ordinates is illustrated in Fig. [l] 



into 



Dvx — 2^loVy = — (c*i + qr d'y + lip] cosTd'z)h, 



T)Vy + (2 — q)^loVx 



-d'h, 



-|- li/^irio sin T = —Qqz' — d'^h, 



(38) 

(39) 
(40) 



Dh : 



-Cs[(9x + qrd'y + lip] cosTd'z)vx + d'yVy + d'^Vz]. (41) 



An alternative form of equation (|4ip . in which the conser- 
vation of mass is manifest, is 



d'tp + d'^ipvx) + dy[p(vy + qrvx)] 

+ d'z[p{Vz + COSTV:^)] = 0. 



(42) 



These equations contain r explicitly but not x' or y' . This 
shows that the local model of a warped disc is horizontally 
homogeneous. As in the standard shearing sheet, every point 
in the x'y' plane is equivalent, if allowance is made for an 
appropriate change of frame of reference. 

We can see from these equations that d = is not 
a solution in the presence of a warp. Equation (|40p would 
require hydrostatic equilibrium in the form d'zh — —Q.'^z , 
but this would provide an unbalanced horizontal force in 
equation (|38|l , as illustrated in Fig. [5] 

The total energy equation associated with equations 

d't{p£) + d'^[(pE +p)v^] +dy[[p£ p){vy + qrv^)] 

+ d'z[{p£ +p){Vz + IV'I COSTVx)] 

= pvx{qO,oVy — \^\Q,Q sin t Vz + IV'l^o cos tQ,oz), (43) 
where 



(44) 
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Figure 5. Regions of high (H) and low (L) pressure in the lo- 
cal representation of a warped disc, showing how the hydrostatic 
vertical pressure gradient leads to a horizontal force that is un- 
balanced by gravity (arrows). The xz plane (of unwarped coordi- 
nates) is depicted at t = 0, with the dotted line corresponding to 
2 = 0. These horizontal pressure gradients oscillate at the orbital 
frequency as the local geometry tilts back and forth, giving rise 
to oscillatory planar motions. 



We will interpret the source term on the right-hand side of 
equation (|43|l below. 



3.4 Fluxes of mass and angular momentum 

In order to connect the local model with the large-scale dy- 
namics of the warped disc described in Section[51 we consider 
the outward radial fluxes of mass and angular momentum. 
The mass flux at the reference radius ro is 



{pUx)hro dr dz' 



(45) 



(in which can be replaced by v^), where the notation {•)h 
denotes horizontal averaging over the coordinates x' and y' 
(where necessary), and the z' integral is over the entire ver- 
tical extent of the disc. The r integral, which in the local 
model is with respect to time, can be interpreted as an inte- 
gral with respect to azimuth as the reference point traverses 
its orbit. Similarly, the angular-momentum flux in the local 
model is (cf. equations 1121 and I13|) 



-roflozea:])^^r()dTdz'. (46) 

As in equation (|12|) . the first two vectors inside the square 
brackets, contributing to the specific angular momentum of 
the fluid, involve non-radial motions combined with the long 
radial lever arm, while the third expression involves the large 
azimuthal motion combined with a vertical lever arm; again, 
we choose not to include the larger constant term involving 
the large azimuthal motion combined with the long radial 
lever arm, which generates the larger flux J-oHq. Writing u 
and z in terms of v and z' , and using the definition (|45|) of 
J-Q, we find that the terms proportional to x on each side of 
the equation cancel, leaving 



(47) 



Go = J J (go)hrodrdz', 
where 

9o ^ pVx{roVy — ro«z By — roQoz' e^) 



is the internal torque per unit area. We then recognize the 
right-hand side of equation (|43[) as the scalar product —[qq ■ 
So)/ro, which is the rate at which energy is extracted, per 
unit volume, by the torque acting on the orbital shear. 

So far we have considered inviscid hydrodynamics. More 
generally, there may be shear stresses in the fluid resulting 
from viscosity or magnetic fields. In the presence of a sym- 
metric stress tensor T, which can also represent the effects 
of turbulence or self-gravitation, the torque per unit area is 
given by the more general expression 

— — pvx{vy fiz — VzSy — Qoz' Bx) — Txy Cz + Txz e.y. (49) 

In the case of a viscous stress, for example, 

Txy = p\-q^ii + + qr 4- jV'l cos r d'^)Vy + d'yVx\ (50) 

and 

Txz = M[|V'|f^o sinr+(a^ + gr + cos r a!,)^^ -f a^f^c] (51) 

are the relevant components, where p is the dynamic viscos- 
ity. 

Since g^-la — goz , the vertical component of the torque 
is given by 



Go ■ Iq 

2^4 



[pVxVy 



Txy) dz' 



(52) 



where the notation (•)h,T denotes horizontal averaging over 
the coordinates x' and y' (where necessary) and averaging 
over the orbital phase r. This is easily understood as involv- 
ing the radial transport of azimuthal momentum (or verti- 
cal angular momentum) by Reynolds and viscous (or other) 
stresses. 

In computing the horizontal components of the torque 
we must bear in mind that the basis vectors ex and de- 
pend on orbital phase. We make use of a complex notation 
and equation (|9]). Thus 



■ (mo 4- ino) = {gox + i5oy)e'^ 
and so 

Go ■ {mo + ino) 



(53) 



27rrg 



[pvx{—^loz' — ivz) + iTxz] dz' 



(54) 



The horizontal components of the torque are conveniently 
encoded in the real and imaginary parts of this expression. 
Not surprisingly, it involves the radial transport of vertical 
momentum (or horizontal angular momentum). 

It is natural to represent the torque for an isothermal 
disc in the form (cf. IOgilvie|[l999l ) 



-27rEcsr (Qi I + QiW + Q'-iW 



]il + Q2r 



dl 



3rl X 



dl 



dr dr 
where Qi, Q2 and Q3 are dimensionless coefficients and 



(55) 



E = 



(56) 



(48) 



is the horizontally averaged surface density of the disc, which 
is independent of r by virtue of mass conservation. This 
representation is just a projection of the torque on to the 
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local basis defined by the geometry of the disc0 Comparison 
of the vertical component of the torque shows that 

-QiEc^ = (/ y" {pv^vy - T^y) dz'^ , (57) 

while comparison of the horizontal components shows that 

^Q4ij\Ec! = (e'^ j [pv^{-noz' - iu,) +ir,,]dz'^ , (58) 

where Q4 — Q2 + iQs is a useful combination. Roughly 
speaking, the Qi term is similar to the usual torque in an 
accretion disc, and (if negative) causes the mass distribution 
to evolve diffusively, while the Q2 term (if positive) causes 
a diffusion of the warp and the Q3 term causes a dispersive 
wavelike propagation of the warp. 

The correspondence between these dimensionles s torque 
coeffic ients and the viscosity coefBcients ui and U2 of lPringlg 
1I1992I) (deriving from the 'naive approach' in Section 2 of 
iPapaloizou fc Pringld[l98^ ) is 

QiCs = -quiD,, Q2cl = \v2^, (59) 

while Q3 has no counterpart in that description. (It would be 
potentially misleading to associate Q3 with a 'viscosity' co- 
efficient because the Qs term is non-dissipative. However, 
the combination Q4 = Q2 + iQa does emerge naturally as a 
complex diffusion coefficient, especially in linear theory.) 

Since the shape of the disc is defined in the local model 
only by the dimensionless warp amplitude | -(/)], it is natu- 
ral to expect the dimensionless torque coefficients Qi, Q2 
and Qz to emerge as functions of as well as any other 
relevant dimensionless parameters such as a. We will see in 
Section |43] how this works for laminar flows, but in the pres- 
ence of instability and turbulence these functions must be 
determined by means of numerical simulations. 



3.5 Computational considerations 

The hydrodynamic equations p8p - (|4H) or their equivalents 
can be solved numerically by standard methods, for example 
by using finite differences on a regular grid in the coordi- 
nates {x',y',z') and imposing periodic boundary conditions 
in x' and y' and some other appropriate boundary condi- 
tions in z' . The additional terms proportional to \tp\ would 
of course require some rewriting of existing codes. Note that 
the primed coordinate grid undergoes an inexorable shear- 
ing, even in the absence of a warp. The coordinates should 
therefore be remapped periodically to avoid excessive distor- 
tion of the grid. This procedure is equivalent to periodically 
resetting the arbitrary origin of time, and can be done with- 
out interpolation if the remapping frequency and the aspect 
ratio of the grid are chosen correctly. In the case of a warped 
disc with its oscillating local geometry, it may be convenient 
to remap once per orbital period, for example by letting r 
run from —n to vr repeatedly, giving rise to a periodic dy- 
namical system. 

In the standard shearing box for an unwarped disc. 



In the absence of a warp, vanishes and m and n are unde- 
fined, but in this case we expect Q to be parallel to I by symmetry. 



there are in fact two different methods of solving the hy- 
drodynamic equations. The more common method is to dis- 
cretize the equations (with time-independent coefficients) in 
a cuboidal domain in unsheared coordinates (x, y, z) and to 
apply (time-dependent) modified periodic boundary condi- 
tions. The other method is to discretize the equations (with 
time-dependent coefficients) in a cuboidal domain in sheared 
coordinates (x', y' , z') (with periodic remapping of the grid) 
and to apply (time-independent) periodic boundary condi- 
tions. For a warped disc, only the second approach is possi- 
ble: the equations must be discretized in primed coordinates, 
because a fixed cuboidal domain in unprimed coordinates 
cannot represent the regions above and below the warped 
midplane of the disc in a way that allows the radial bound- 
aries to be identified (cf. Fig. [3). This limitation is related to 
the fact that vertical gravity is essential to the description 
of a warped disc, and the z direction cannot be regarded as 
periodic. When supplied with periodic boundary conditions 
in x' and y' , our local model takes the form of a warped 
shearing box. 

4 LAMINAR FLOWS 
4.1 Separation of variables 

The simplest solutions of equations (|38|) - (|4ip are in fact 
independent of x' and y' and 27r-periodic in r: they are hor- 
izontally uniform and oscillate at the orbital frequency. (In 
a global context, such local solutions correspond to the sit- 
uation in which the warped disc has a structure that varies 
on a horizontal length-scale much larger than the thickness 
of the disc, and on a time-scale much longer than the orbital 



time-scale.) These laminar flows satisfy 

Dvx — 2^loVy = —\^\cosT dih, (60) 

T)vy + (2 - q)noV:c = 0, (61) 

T)Vz + l^plflo sinr Vx = —i^oz' — d'zh, (62) 

T>h = —Cs{\^\cosT d'^Vx + d'zVz), (63) 

with 

D = d't + {v, + \^p\cosTVx)d',. (64) 
A nonlinear separation of variables is then possible, with 

Vx{z',t') = u{t)D,oz', (65) 

Vy{z' ,t') — v{t)Qoz' , (66) 

Vz{z' ,t') — w{t)Qoz' , (67) 

h{z',t')^c!f{T)-^Qlz'^g{T), (68) 
where u, v, w, f and g are all dimensionless. Thus 

drU + (to + \tp\ cos T u)u — 2u = j?/;] cos r g, (69) 

d-rV + (to + \^l>\ cos r u)v -|- (2 — q)u — 0, (70) 

d^TO + (w + \ tp\ cos r u)w + \tp\ sin r u — g — 1, (71) 

d^/ = — (to + I?/)] cos r It), (72) 

d^p = — 2(to -I- cosrit)(7, (73) 
where d^ denotes the ordinary derivative d/dr. 
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When a dynamic shear viscosity fi = ap/Q,o and a dy- 
namic bulk viscosity ~ ahp/flo are included, where a and 
Ob are constant dimensionless coefficients (thus providing a 
kinematic shear viscosity u = aCs/Oo, etc.), these equations 
become 



drtt + (w + 1-01 COS T u)u — 2w = j-j/jj COS r g 

— (ah + cos T g{w + \xp \ cos r u) 
— a(7[|?/>| sinr + (1 + It/j]^ cos^ 

dTi" + (to + I?/"! cos r it)i; + (2 — g)it 

= —ag[—q\'4)\ cost + (1 + l?/"!^ cos^ 

drW + (w + I V'j cos r u)w + I^/jI sin t u = g — 1 

— (Qb + |a)(7(TO + IV"! cos r it) 

— a(7[|?/>|^ sin r cos r + (1 + cos^ 1")^], 



drf = —{'w + \tp \ cos r u), 
drfl = — 2(to + cos r it)(7. 



(74) 
(75) 

(76) 

(77) 
(78) 



Our immerical treatment of these ordinary differential 
equations in Section [4.31 below provides us with a family of 
'exact' nonlinear solutions of the local model representing 
the laminar flows in a warped disc. The same laminar flows, 
combined with the warping motion, can alternatively be re- 
garded as nonlinear solutions of equations (|21|| - (|24p in the 
variables of the standard shearing sheet, although they do 
not then appear to be horizontally uniform and it would not 
be obvious how to obtain them without making the trans- 
formation to warped coordinates. For example, the vertical 
velocity would be 



Uz ~ lipl^io sinr a; + ui(r)f2o(^ + |^| cos r x). 



(79) 



The differential equations for / and g are related in 
such a way that the surface density of the disc, which is 
proportional to g~^^^, is independent of r. The total en- 
ergy equation for the laminar flows has the form 



d. 



(u^ + + + p In g -|- 1) 



- —{quv — \tp\ sinr MUi + cos r it) 

-a]?/'! (sin T u — q cos t v -'r\il)\ sin r cos r w) 
-a(l -|- jT/il^ cos^ + v'^ w^) 

-{at, + ha){w + cos tu)'^ . 



(80) 



The first expression on the right-hand side is equivalent to 
the right-hand side of equation (|43|l and corresponds to the 
extraction of energy by the torque acting on the orbital 
shear. The third and fourth expressions are negative definite 
and correspond to viscous damping of the fiow. The second 
expression is proportional to viscosity but is not negative 
definite; its interpretation is not straightforward. 

In Appendix[X]we give the corresponding equations for 
ad iabatic flow. T hey are exactly equivalent to those derived 
bv lOgilvid l| 19991 ) in a different way, by writing the hydro- 
dynamic equations in a warped spherical polar coordinate 
system that follows the global distortion of the warped mid- 
plane, and carrying out an asymptotic expansion of the so- 
lution in the limit of a thin disc. In that paper the orbital 



phase r is directly related to the azimuthal angle 4>. Here 
we have not explicitly made use of an asymptotic expan- 
sion, but have derived the equations consistently within the 
context of the local approximation. 



4.2 Horizontal and vertical oscillators 

The physical interpretation of the laminar flows is assisted 
by changing to the variables 

{U,V,W) ^ {u,v,w)H, H^g-^''^. (81) 
Since the density and pressure then have the form 



/5 oc p oc exp < /(r 



(82) 



it can be seen that H{t) is the scale- height of the disc in 
units of the scale- height Cs/ilo of an equilibrium unwarped 
disc. The variables (U, V, W) are therefore more like mo- 
menta, or vertically integrated velocities. 
In the inviscid case we then have 



(83) 
(84) 
(85) 
(86) 



d^f/-2y = lVlcosr//"\ 

drV + {2 ~q)U = 0, 

d^W + \ip\ sinr 17 = R-^ - H, 

d^H = W +\il>\cosrU. 

In the absence of a warp {\ip\ =0), equations (I83|l and (|84|l 
describe a linear epicyclic oscillator. They can be combined 
to give d^V — — 2(2— g) V, or the same equation for U. Equa- 
tions (|85|) and (|86|) describe a nonlinear vertical oscillator, 
the 'breathing mode' of the disc. They can be combined 
to give drH — — H. The vertical oscillator therefore 

has a potential energy function — IniJ" -I- + constant, 
representing the sum of internal energy and gravitational 
energy, with equilibrium at _ff = 1. It can in principle sup- 
port oscillations of any energy, but this may involve severe 
compression because of the slow divergence of the potential 
as -ff — >■ 0. The natural frequency of the oscillator in the 
linear regime is \/2 (or -^7 -f 1 for an arbitrary adiabatic 
exponent 7); more generally, it depends on the amplitude of 
the oscillation. 

The horizontal and vertical oscillators are coupled in 
the presence of a warp. The |i/'| term in equation (|83|l rep- 
resents the radial pressure gradient that occurs in a warped 
disc: in the absence of a warp, the pressure varies with 2 
to provide hydrostatic equilibrium, and warping of the disc 
causes a radial gradient to appear (Fig. The |?/>| term 
in equation (|85p represents the horizontal transport of the 
vertical momentum of the warped orbital motion. Finally, 
the \tp\ term in equation (|86p represents a contribution to 
the velocity divergence from horizontal motion in a warped 
disc. These coupled horizontal and vertical oscillators fea- 
tured prominently in the analys is of weakly no nlinear bend- 
ing waves in Keplerian discs bv lOgilvid (|2006l ). 

As we have seen, the viscous terms in a warped disc have 
a complicated form. In the absence of a warp, the equations 
for laminar viscous flows can be written 



drU -2V = -aH'^U, 
drV +{2-q)U = -aH-'^V, 



(87) 
(88) 
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drW = H-^ -H-{ab + ^^a)H-^W, 
drH = W, 



(89) 
(90) 



showing that the horizontal motion is damped by shear vis- 
cosity while the vertical motion is damped by both shear 
and bulk viscosity. 

4.3 Numerical solutions 

Typical numerical solutions for the laminar flows are shown 
in Figs |6] and [T] These were computed by a standard shoot- 
ing method using a Runge-Kutta integrator with adaptive 
stepsize. We first consider the inviscid case a = Ob = for 
a non-Keplerian disc with g = 1.6 (Fig.|6] panels 1-3). Note 
that the regime | < q < 2, in which the epicyclic frequency 
is positive but less than the orbital frequency, corresponds 
to conditions close to a black hole, where relativistic effects 
are important. For sufficiently small warp amplitude the 
horizontal flows u and v are proportional to and have 
a sinusoidal dependence on orbital phase, corresponding to 
the linear horizontal oscillator being driven at the orbital 
frequency, which here is greater than the natural (epicyclic) 
frequency of the oscillator (panels 1-2). The vertical veloc- 
ity w and the departure from hydrostatic equilibrium, g — 1, 
are proportional to and therefore smaller. (Note that 
w does not include the vertical velocity associated with the 
warped orbital motion itself.) 

As the warp amplitude is increased further, however, 
the behaviour of the solutions changes markedly and the 
branch of inviscid solutions terminates; this happens at 
\il>\ ~ 0.261 in the ca se q = l . Q (pan el 3). This phenomenon 
was already noted bylOgilyj^ l| 19991 ') and is shown in fig. 2 of 
that paper. In Appendix iBl we explore mathematically the 
breakdown of the inviscid solution and find that this occurs 
at « (<J - |)/l-45 for < (g - |) < 1. The physical 
reason for the breakdown is a type of nonlinear resonance 
involving the coupled horizontal and vertical oscillators. If 
a viscosity is included, it is possible to avoid the breakdown 
and obtain solutions for larger values of but they can 
have quite extreme properties (panels 4-6). 

For non-Keplerian discs with g < | no such breakdown 
occurs and the inviscid solutions can be continued to large 
values of It/; I (Fig.[7l panels 1-3). The phases of the horizon- 
tal velocities differ by tt from the case g < § because the 
orbital frequency is now less than the epicyclic frequency. 

The Keplerian case q = | is of course of greatest inter- 
est. Here no 27r-periodic solutions can be obtained without 
introducing a viscosity to moderate the resonance resulting 
from the coincidence of orbital and epicyclic frequencies. For 
sufficiently small the horizontal velocities u and v are 
sinusoidal and proportional to I'i/'l/c*; their phases are in- 
termediate between those of the inviscid fiows in the cases 
g > I and g < I because the horizontal oscillator is being 
driven at its natural frequency, rather than above or below 
it (panels 4-5). The solutions can be continued to large val- 
ues of IV' I, but again they can have quite extreme properties 
(panel 6). 



4.4 Role of time-dependence of the warp 

In constructing the local model of a warped disc we have 
assumed that the geometry of the warp is fixed in a non- 



rotating frame of reference. If, however, the warp takes the 
form of a propagating bending wave or a rigidly precess- 
ing structure, for example, then this time-dependence may 
have an important effect on the local dynamics even if the 
time-scale for the evolution of the warp is long compared to 
the orbital time-scale. This is because the horizontal forc- 
ing due to the warp no longer occurs at exactly the epicyclic 
frequency in a Keplerian disc, and the exact resonance is bro- 
ken. The easiest way to model this effect within the present 
framework is to allow the parameter q to be adjusted from 
its resonant value of |. If, for example, the warp precesses 
slowly at a rate fip (< for retrograde precession), which is 
determined by global considerations, then in the fiuid frame 
the geometry oscillates at frequency 0.-0.^. If the disc is 
Keplerian, then the required detuning of fip between the 
horizontal forcing frequency and the epicyclic frequency can 
be achieved by adjusting q such that [\/2(2 — g) — l]f2 = f2p. 

4.5 Torques 

In Section lH^ we showed how to compute the internal torque 
in the local model. For the laminar viscous flows described 
in Section [4. II the dimensionless torque coefficients evaluate 
to 



Qi = {—g ^uv + a{—q+ lit] COST v))t 
and 



(91) 



(54|'!/'| — {e'^[g "'"ii(l-l-iw)— iad?/;! sinr-t-ji/jj COST to-|-m)])t.(92) 

T hese are the di rect analogues of equations (112) and (120) 
in IOgilvi3 (jl999l ). but are expressed in a different notation, 
and are slightly simplified because the disc is isothermal. 
The case of adiabatic fiow is treated in Appendix 1X1 

Note that, if the laminar fiow is neglected by setting 
u = V = w = 0, we obtain the simple results Qi = — aq 
and Q4 = ia. These torque coefficients represent only the 
viscous stresses associated with the warped orbital motion, 
and would lead to a diffusion of the warp on (twice) the 
viscous timescale as found erroneously in early theoretical 
work on warped discs. In practice the torque (mainly from 
the radial advective transport of horizontal angular momen- 
tum) associated with the laminar fiows exceeds the viscous 
torque, and leads to a more rapid evolution of the warp. 

For an inviscid disc, which must be non-Keplerian for 
laminar solutions to exist, only Q3 is non-zero and it deter- 
mines the (dispersive) propagation of bending waves accord- 
ing to the equation 



(93) 



which follows from equation Q with 7^ = in this case. 
The variation of Qs with and q is sho wn in Fig. O which 
is equivalent to fig. 2 of lOgilvid (|l999t ) except that it is 
for an isothermal disc. As noted in Section [4.31 the inviscid 
solutions break down for | < q < 2 when \ip\ exceeds a 
critical value, owing to a type of nonlinear resonance. The 
analytical result for sufficiently small \tp\'^/{2q — 3) is 



1 



2(2q - 3) 



1 + 



+ 



(94) 



(2q-3) \{2q-3)\ 

In Fig. [9] we show the variation of all three coefficients 
with |-i/)| and a for a Keplerian disc without bulk viscosity. 
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Figure 6. Selection of laminar flows for g = 1.6. Solid (black), dotted (red), dashed (green) and dot-dashed (blue) lines represent u, v, 
w and g — 1, respectively (see equations I65lf68t . The inviscid solutions terminate at [ipl ~ 0.261 in this case, but larger values of \ip\ can 
be reached by including viscosity. 
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Figure 7. Continuation of Fig. [6]for q = 1.4 and q = 1.5. The inviscid solutions for q = 1.4 continue to large values of \^\. For q = 1.5 
some viscosity is required. 
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Figure 9. Dependence of the torque coefficients Qi (top), Q2 (middle) and Q3 (bottom) of laminar flows on the warp amplitude and 
viscosity for a Keplerian disc without bulk viscosity. The right-hand panels enlarge the region < a < 0.2. 
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Figure 8. Dependence of the torque coefficient Q3 of laminar 
flows on the warp amplitude and dimensionless shear rate for an 
inviscid non-Keplerian disc. In the blank sector (top right) no 
solutions are found. The short unlabelled contours have values 
±10. 



These results are equivalent to figs 3-5 of lOgilvid ()l999l ) 
except that they are for an isothermal disc. The analytical 
results for sufficiently small in this case are 



Oi 



and 



3a 
2 



(1 - 8q^ + 3a^) 
4Q(4 + a2) 



V'|^ + o(lv:.|'*) 



Q 



1 + 2ia + 



2a{2 + ia) 
which can be decomposed into 
l + Ta^ 



Q 



a(4 + a^] 
3-6q^ 



-om 



(95) 



(96) 



(97) 



'53^2(4f^+^(l'^>')- (98) 

The figures show that the variation with amplitude is mild 
except when the viscosity is small, as highlighted in the 
right-hand panels. In particular, the coefficient Qi reverses 
sign for sufficiently large and small a. This angular- 
momentum flux reversal occurs when the Reynolds stress 
associated with the correlation of u and v exceeds the usual 
viscous stress. If it occurs, it can be expected to cause an 
antidiffusion of the mass distribution of the disc, leading 
to the formation of disjoint rings0 The warp diffusion co- 
efficient Q2 diverges at small a and small \'4>\ where the 
resonance is located, but this effect is weakened by nonlin- 
earity at larger amplitude. As shown in Paper II, however, 
the laminar solutions are generally unstable at small a and 

Angular-momentum flux reversal occurs in planetary rings 
that are strongly perturbed by nonlinear density waves, and 
is thought to be responsible for the f ormation of sharp edges 
teorderies. Goldreich fc Tremainelll989l) . 



large so the torque coefficients cannot be trusted in this 
region. 



5 CONCLUSION 

This paper offers several advances in the theory of warped 
discs. First, we have const ructed a local model that genera l- 
izes the shearing sheet of iGoldreich fc Lvnden-Belll (|l965h . 
The warped shearing sheet is horizontally homogeneous and 
admits periodic boundary conditions, although its geometry 
oscillates at the orbital frequency. This leads to a computa- 
tional model in the form of a warped shearing box, which 
can be used for linear and nonlinear studies of fluid dynam- 
ics, magnetohydrodynamics, etc. in warped discs, especially 
to investigate processes such as instability and turbulence 
that cannot be resolved in global numerical simulations. It 
is related to the standard shearing box by a relatively simple 
coordinate transformation described in Section [3.31 Second, 
we have shown how to use the local model to calculate all 
three components of the torque that governs the large-scale 
evolution of the shape and mass distribution of the disc. 
Third, we have provided an independent and simpl er route 
to the nonlinear theory of warped discs derived by IOgilvi3 
(|l999l ), showing how the solutions obtained in that paper 
correspond to the simplest laminar flows that can occur in 
the local model. In Paper II we use the local model to anal- 
yse the linear hydrodynamic stability of these laminar flows. 
The widespread instability that we flnd there requires fur- 
ther investigation in future work. 
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APPENDIX A: ADIABATIC FLOW 

The basic equations governing the adiabatic flow of an in- 
viscid fluid in the local model are 



2 1 

Dux — 2f2oWy = "igriox dxP, 



--dyP, 
P 



2 1 

Dm^ = —Q.QZ dzP, 



Dp = -p{dxUx + dyUy + dzUz), 

Dp = -yp{dxUx + dyUy + dzUz), 
where 

D = 9t + Uxdx + Uydy + Uzdz 



(Al) 
(A2) 

(A3) 

(A4) 
(A5) 

(A6) 



is the Lagrangian derivative and 7 is the adiabatic expo- 
nent, which we assume to be constant. By introducing the 
speciflc internal energy e, the specific enthalpy /i, the specific 
entropy s and the temperature T, which satisfy the relations 



h = e + 



P 



Ah = Tds + 



dp 



the basic equations can be rewritten in the form 

Hux — 2Q.oUy = 2qQ^x — dxh + TdxS, 

T)Uy + 2QoUx = —dyh + TdyS, 

Duz = —Clo^ — dzh -\- TdzS, 

Bh = —(7 — l)h{dxUx + dyUy -h dzUz), 

Ds = 0, 



(A7) 
(A8) 

(A9) 
(AlO) 
(All) 
(A12) 
(A13) 



assuming an ideal gas for which h = 'ye. In warped shearing 
coordinates, these equations take the form 

T)Vx — 2QoVy = —{d'x + qrd'y + \'ip\ cost d'z)h 

+T{d'x +qTd'y + \ip\ cosTd'z)s, (A14) 



Bvy + {2- q)noVx = -dyh + Td'yS, (A15) 
Bvz + \^\no sinrvx = -0.lz' - d'zh + Td'zS, (A16) 

Bh = -(7 - l)h[{d'x + qrd'y + \tp\ cosTd'z)vx + d'yVy + d'zVz], 

(A17) 



Ds = 0, (A18) 
with 

B = d't + Vxd'x + {vy + qTVx)d'y + {vz + cos r Wa:)9^,(A19) 

where v is the relative velocity. Laminar fiows independent 
of x' and y' satisfy 

Bvx -2noVy = -\ip\ cos T d'zh + T\il;\ cos rd'zS, (A20) 
Bvy + {2- q)(loVx = 0, (A21) 
Dt;;^ -I- IV'ino sinrVx = -fio^' - d'zh + Td'zS, (A22) 
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Dh = -{j - l)h{\iP\ COST div^ + d',v,), (A23) 

Ds = 0, (A24) 
with 

D = a; + (w. + iV'|cosru,)a;. (A25) 

A nonlinear separation of variables is then possible, with 

v^{z',t') = u{r)noz', (A26) 

Vy{z' ,t') ^ v{t)Qoz' , (A27) 

v4z',t') =w{T)noz', (A28) 

hiz',t')^nl[fir)-^gir)z% (A29) 

s{z',t') ^ s{t), (A30) 

where u, v, w and g are all dimensionless. Thus 

drU + (w + \tp\ COS T u)u — 2v = \'4>\ COS r g, (A31) 

d-rV + {w + \'4)\ COST u)v + {2 - q)u = 0, (A32) 

dr^y + (™ + IV'I COS r 11)10 + I'f/'l sin t u = g — 1, (A33) 

d^/ = -(7-l)(TO + lVlcosru)/, (A34) 

dy-g — —{j + l){w + \ip\cosT u)g, (A35) 

drs = 0. (A36) 

When viscosity is included in the equation of motion, but 
viscous heating is neglected, these equations become 

drU + (w + \ip\ COS T u)u — 2v = \tp\ cos T g 
— (ob + |a)|V'l cos T g{w + \tp\ cos t u) 
—aglltpl sinr + (1 + cos^ 



where 



drV + (to + \tp\ cos T u)v + (2 — q)u 

= —ag[—q\'4)\ cost + (1 + l?/"!^ cos^ 



(A37) 
(A38) 



drW -\- {w -\- cos r u)'w + I?/)! sin r it = g — 1 
— (ab + + IV;! cos r ii) 

— Q(?[|?/'|^ sin r cos r + (1 + j-!/"!^ cos^ t)w\, (A39) 



drf = — 1)(™ + l?/^! COS T u)f, 

drg = —(7 + l)(if + IV"! cosTu)g, 

drS = 0. 



(A40) 
(A41) 
(A42) 



These are ex actly equivalent to equations (105)-(109) of 

iQgiivid (Ulii). The mapping from the earlier notation to 

the current notation is as follows: /i 1— > /, /2 i-> g, fa u, 
fi ^ —{w + \'ip\ costm), /s 1-^ i;, r h> 7, ft^ i-^ 2(2 — g), 
(j) T. The equations for / and g are related in such a way 
that the surface density of the disc, which is proportional to 
^(7+i)/2(7-i)g-i/2^ is independent of r. 

The dimensionless torque coefficients for a non- 
isothermal disc must be defined differently fro m equa- 
tion (1551) . because Cs is no longer constant. As in lOgilvid 
l|l999l ). we define instead 



g = --2'Kir'^n^ [Qii 



. dl ^ , dl 
)2r-- + Qarl x 



dr 



dr 



(A44) 



is the orbitally averaged second vertical moment of the den- 
sity. Prior to averaging, its variation with orbital phase is 
proportional to oc fg~^, and we define /g — fg~^/{fg~^)r 
as in the earlier notation. Then the dimensionless torque 
coefficients for these laminar flows are given by 



= {fe[-uv + ag{-q+ costv)])^, 



(A45) 



Q4,\tp\ = (e'^/6[u(l + iw) — iag{\tp\ sin t + \^()\ cos tw + u)])t, 

(A46) 

which agree with equations (112) and (120) in lOgilvid 
(,1999,1 . 



APPENDIX B: EXISTENCE OF INVISCID 
LAMINAR FLOW SOLUTIONS 



Equations (|83p - (|86|) for inviscid laminar flows are 

d^t/-2y = lVlcosr_f/'"\ (Bl) 

d^V + {2 - q)U = 0, (B2) 

d^W^+|7/>|sinrf/ = //"^-_H", (B3) 

d^H = W +\iP\costU. (B4) 

Note that U and W are expected to be odd in r while V 
and H are even. 

The regular expansion of the solution for small \tp\ is 

U=\i,\Ui+0{W'), (B5) 

V ^WVi + 0{\ipf), (B6) 

W = \i^fW2 + 0{\i)f), (B7) 

H = 1 + W^H2 + 0{W''). (B8) 

At leading order we find 

d.,;7i - 2Vi = COST, (B9) 

d,Vi + (2-g)(7i =0, (BIO) 
with the solution (having the required period of 2n) 
sin T 



Ui = 
Vi = 



2g-3' 

(2 — q) cos T 
2q~3 ' 



(Bll) 



(B12) 



provided that g 7^ |. This solution breaks down in the Ke- 
plerian case because of the coincidence between the orbital 
frequency and the epicyclic frequency. 

Solutions valid for small \ip\ and q close to | can be 



(A43) 



found instead by the following expansion: 

U ^\ij\'^[Uo + WU2 + 0{\iP\^)], 

V ^W-\Vo + \i>\^V2 + o{W% 

W = Wo + 0{\^\^), 



(B13) 
(B14) 
(B15) 
(B16) 
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Then we find 
d,f/o- 21/0 = 0, 

drVo + ic/o = 0, 



d^Wo + sin tUo = Ho^ - Ho, 

drHo = Wo+ COS T Uo, 
drf/2-2V2 = COST Ho \ 



d^V2 + -U2 - QUo = 0. 



(B17) 

(B18) 
(B19) 

(B20) 
(B21) 
(B22) 

(B23) 



Equations (fBT8|) and IfOW} give (d? + l)((7o, Vb) = (a 
free linear epicyclic oscillator) and liave the general solution 



0.8 



0.6 



0.4 



0.2 



0.0 




Vo = A cos T + B sin r, 
C/o = 2A sin r — 2B cos r. 

We set B = on grounds of symmetry. 
Equations (|B20|l and (IB21|l give 



Ho 



Ho + A(3cos2r - 1). 



(B24) 
(B25) 



(B26) 



Here we see the nonlinear vertical oscillator forced by the 
coupling (through the warp) to the horizontal oscillator. 
The natural frequency of the vertical oscillator in the linear 
regime is \/2, which is far from resonance with the forcing 
at frequency 2. However, the natural frequency depends on 
the amplitude of the oscillation. The constant term —A in 
the forcing also shifts the centre of the oscillation away from 
the equilibrium position Ho ~ 1. 

Equations (lB22l) and (|B23ll give 



(d^ + 1) V2 = 2QA cos r - - cos t Ho\ 



(B27) 



a forced linear epicyclic oscillator. For a periodic solution to 
exist, the right-hand side must contain no Fourier compo- 
nent proportional to cost. The mean value and the cos2r 
content of H^^ depend nonlinear ly on A. The aforemen- 
tioned solvability condition is therefore a nonlinear equation 
for A, which may or may not have a solution. 

For sufficiently small A, the solution of equation (|B26|) 
can be obtained by a series expansion. It is 



Ho 



and so 



-^(l + 3cos2r)-f^(77- 
+0{A'), 



I cos 2r — 9 cos 4r) 
(B28) 



2QA cos ~ 2 ''^^ ^ ^0 ^ ~ ~ 2 '^^^ ^ 

-f — [(16(5 ~ 5) COST — 3cos3t] 
8 
42 

(406 cos r + 387 cos 3r -f 135 cos 5r) 

448 

-fO(A'). 



Figure Bl. The numerically determined function of the ampli- 
tude A that should equal QA when the solvability condition for 
equation I IB27II is applied. The black solid line, which terminates 
at A ^ 0.36, shows the component of cos r in the Fourier series of 
(1/4) cos tH^^ determined from the numerical solution of equa- 
tion l)B26[l . The red dashed line shows the series approximation 
given by the right-hand side of equation I IB30II . The blue dotted 
line is of slope 1.45 and passes through the origin. 



This series can be continued to higher order and is accurate 
if \A\ is not too large. The solvability condition gives 



«^=4 + T6 + 



29A^ 



+ ■ 



333A-^ 



+ 



270913^'' 



64 448 200704 
15848541yl'5 36001767715^"^ 



+ 



5970944 6496387072 
253686124482669^'' , .g. 



20969525420032 



+ 0{A'' 



(B30) 



(B29) 



For larger \A\, we solve equation (|B26|) numerically instead. 
The solution breaks down for A > 0.36, although it can be 
continued to large negative A. (In fact the solution of equa- 
tion lB26l is not unique, but we consider here only the solution 
that is consistent with the expansion 16281 ) In Fig. IBll we 
plot the numerically determined function of A that should 
equal QA when the solvability condition is applied. The red 
dashed line shows the series approximation given by the 
right-hand side of equation HB30|I . The solvability condition 
is satisfied when the black line in Fig. IBll intersects with a 
straight line of slope Q through the origin. For negative val- 
ues of Q this condition can always be met for some negative 
value of A, but for positive values of Q a solution exists only 
if Q is sufficiently large. The critical slope is indicated by 
the blue dotted line and corresponds to Q ~ 1.45. Therefore 
there are no solutions for < Q < 1.45. In the full numerical 
problem for inviscid laminar flows, the solution does indeed 
break down at Q « 1.45 and A ~ 0.30. 

The breakdown of the solution can be understood as a 
type of nonlinear resonance. It occurs while the nonlinear 
vertical oscillator exhibits a finite response, but the anhar- 
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monic and asymmetric character of that oscillator is essen- 
tial to the behaviour. 
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